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THE INTERNATIONAL GEOMAGNETIC REFERENCE FIELD I965.O 
' ' ' IN DIPOLE COORDINATES * 

Gilbert D. Mead 
Laboratory for Space Physics 
NASA Goddard Space Flight Center 
Greenbelt, Md. 20771 

At the working group session of the meeting on Quantitative Magnetospheric 
Models held in Boulder March 18-20, 1970, the need was recognized for a ref- 
erence model of the earth’s internal field with input and output in geomagnetic 
dipole coordinates ( Chapman, 1963) instead of geographic coordinates. Such a 
model could be used much more readily in conjunction with various models of 
external field sources (e.g., ring current, magnetopause, or tail fields) 
whose orientation depends upon the direction of the magnetic dipole axis 
rather than the rotation axis. 'Model field and field-line -tracing calculations 
could then be performed in a physically more meaningful coordinate system whose 
z-axis is aligned with the magnetic dipole; simple longitudinal transformations 
would suffice to link this system with the solar magnetic coordinate system, 
which is also aligned with the dipole. Deviations from a dipole field would 
be much more readily identifiable. 

A program has been written to transform the spherical harmonic coefficients 
of any model field into any arbitrarily-tilted coordinate system, using a method 
similar to Stera^s (1965). For each value of n in the spherical harmonic 
expansion, the scalar potential V was calculated at 2 n + 1 randomly-selected 
points and set equal to the potential in the tilted system. The resulting 
equations were solved for the 2 n + 1 unknown coefficients. The process was 
then repeated for the time derivatives of the coefficients representing the 
secular change. 

* To be published (without Fortran listings) in the Journal of Geophysical 
Research , Vol. 75 > 1970* 
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The resulting coefficients for the International Geomagnetic Reference 

Field I 965 .O ( IAGA Commission 2 Working Group 4, I 969 ) are given in Table 1. 

The normalisation of the Legendre polynomials and all definitions are the 
IAGA 

same as in the . reference except for the orientation of the coordinate system. 

/\ . 

1 

The geographic coordinates of the new pole are: 

0 Q = 11.455° geocentric colatitude 

\ o = - 69 . 76 I 0 east longitude 

as determined from the three n = 1 coefficients of the IGRF field at 
epoch 1965 .O. The input to a subroutine using these coefficients is 
geomagnetic dipole colatitude ( 8^) with respect to this pole and east 
geomagnetic dipole longitude (X^) with respect to a meridian passing through 
the north and south dipole poles and the south geographic pole. The 
(northward) and Y^ (eastward) components of the field intensity obtained 
as output are then measured with respect to the north dipole pole rather than 
the geographic pole. The Z (downward) component is independent of the 
orientation of the coordinate system. 

It is often necessary to transform positions and/or field components 
from one system to the other. The geometry is shown in Figure 1. P is 
the geographic pole, D is the dipole pole, G is Greenwich, and A is a •, 

point on the reference sphere (or the intersection of a vector through jV 

< 

upon the celestial sphere) at which the field is to be evaluated; 9 and * 

6 

X sire the geographic colatitude and longitude of A; and 6 is the rotational 

D 

i 

angle at A used to transform the X and Y field components; 6 is positive 
if D appears to be west of P, negative if east. In a pure dipole field, 

6 would be the negative of the declination angle everywhere. 


t 


3 



Figure 1. Earth's surface (or, alternatively, surface of celestial f 

sphere), showing spherical triangle DPA needed for coordinate transformations. 

The dipole tilt is exaggerated for greater clarity. I 

| 
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From the solution of the spherical triangle DPA, 

cos0, = cos0 cos0 + sin0 sin0 cos(X -X ) 
d o g o g g o 

sin0^ » + y 1 - cos 2 0^' 

cosX, = (cos0 cos0, - cos0 ) / ( sin0 sin0,) 
a o a 6 o a 

sinX, = sin0 sin(X -X ) / sinQ, 
a g g o a 

cos6 = (cos0 - cos©, cos0 ) / (sin0. sin0 ) 
o d g a g 

sin6 = sin0 sin(X -X ) / sin0. 

o go' a 

and 

X. = X cos6 - Y sin6 

d g g 

Y, a x„ sin5 + Y cosfi 

d g g 

where the sine and cosine of each angle determine the proper quadrant. 

The inverse transformations are given by- 

cos© a cos© cos0. - sin0 sin0, cosX. 
g o d odd 

cos(X -X ) * (cos0 - cos6 cos© ) / (sin0 sin0 ) 

o ^ ^ ^ o ^ o 

sin(X -Xj = sin0, sinX. / sin0„ 
go d d g 

with similar expressions for the rotation of the and Y^ field components. 

The dipole latitude will always differ from geographic latitude by less 

O i 

than 11.5 . Except near the poles, dipole longitude will generally be ) 

about 70 ° greater than geographic longitude and the angle 6 will usually 
be less than ± 15 °« In a pure dipole field, Y^ « 0 everywhere, and in the 
IGRF field, Y^ is small and approaches zero rapidly at large distances. 


5 


r 


The coefficients were checked by calculating the field components 
at a number of points and at several time periods, both in geographic 
coordinates, using the original coefficients, and in the corresponding 
dipole coordinates, using the present coefficients. The calculated 
intensities agreed with each other to within about 4 gammas at the 
earth's surface and 0.2 gammas at a geocentric distance of 3 Rg- 
This is consistent with the fact that each of the transformed coeffi- 
cients was rounded to the nearest gamma (derivatives to the nearest 
tenth gamma/yr). Field directions corresponded to within 0.005 degrees. 

A Fortran program to compute field values in geomagnetic dipole 
coordinates from the IGRF 1965-0 field, as well as a program to transform 
positions and field components from geographic to dipole coordinates and 
vice versa, can be obtained from the author or from the U.S. National 
Space Science Data Center, NASA Goddard Space Flight Center, Greenbelt, 
Md. 20771. 
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Table 1. IGRF I965.O Dipole Coefficients: 
e Q = 11.435°, X Q = -69.761° 


n 

m 

m 

g n 

h m 

n 

• m 
6 n 

h ffi 

n 

1 

0 

-30953 v 


16.0 y/yr 

1 

1 

0 

0 

2.0 

7.4 

2 

0 

-618 


-18.8 


2 

1 

2997 

2255 

20.9 

-1.4 

2 

2 

-1875 

481 

8.8 

12.2 

3 

0 

906 


-3.4 


3 

1 

-1238 

-1758 

-6.8 

-8.0 

3 

2 

-1052 

1170 

1.3 

4.6 

3 

3 

-546 

-485 

-0.5 

8.0 

k 

0 

837 


-0.3 


u . 

1 

-496 

962 

l.l 

-1.3 

k 

2 

15 

176 

1.2 

-3.5 

k 

3 

311 

-39 

-0.3 

-0.7 

k 

4 

-317 

-312 

-4.7 

1.9 

5 

0 

-l4o 


-0.5 


5 

1 

90 

344 

-3.9 

2.1 

5 

2 

-322 

-49 

-2.0 

0.3 

5 

3 

53 

170 

0.2 

1.4 

5 

4 

-138 

103 

1.4 

-0.7 

5 

5 

7 

46 

0.8 

-0.4 

6 

0 

48 


0.3 


6 

1 

-17 

7 

0.0 

-0.4 

6 

2 

60 

-55 

-1.3 

-0.3 

6 

3 

185 

78 

-0.7 

-2.3 

6 

4 

-151 

-56 

-0.9 

1.7 

6 

5 

22 

-29 

-0.2 

-0.2 

6 

6 

-48 

-95 

-0.8 

0.3 

7 

0 

69 


0.6 


7 

1 

-4o 

-39 

1.0 

-1.0 

7 

2 

8 

61 

0.0 

-0.1 

7 

3 

-14 

-9 

0.5 

0.1 

7 

4 

15 

36 

-0.2 

0.0 

7 

1 

5 

3 

4 

-0.1 

0.3 

7 

6 

27 

-3 

-0.2 

-0.5 

7 

7 

5 

10 

0.2 

-0.5 

8 

0 

10 

• 

-0.1 


8 

1 

-1 

12 

-0.3 

0.7 

8 

2 

12 

-2 

-0.3 

0.1 

8 

3 

-11 

1 

0.0 

-0.1 

8 

4 

-14 

4 

-0.1 

-0.2 

8 

5 

4 

7 

0.3 

-0.1 

8 

6 

-26 

10 

0.6 

0.2 

8 

7 

4 

11 

0.2 

0.1 

8 

8 

-10 

9 

0.2 

0.4 
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SUBROUTINE TS FORM ( GLAT, G LO N, D LA T,D LON, SIND, COSD, I ) 

PROGRAM TO TRANSFORM GEOGRAPHIC COORDINATES TO GEOMAGNETIC DIPOLE 
COORDINATES (1*1), OR VICE VERSA (1-2). SEE MEA0(1970). 

INPUT (OR, IF 1-2, OUTPUT) ARE GEOCENTRIC LATITUDE AND EAST 
LONGITUDE IN DEGREES. OUTPUT (OR, IF 1-2, INPUT) ARE DIPOLE 
LATITUDE AND LONGITUDE W/R/T THE DIPOLE MERIDIAN PASSING THROUGH 
THE SUUTH GEOGRAPHIC POLE. SIND AND COSD ARE USED TO ROTATE THE 
X AND Y FIELD COMPONENTS FROM GEOGRAPHIC TO GEOMAGNETIC OR V/V. 

IF (RAO. EQ. 57. 29578) GO TO (1,2), I 

RAD » 57.29578 

COLAT - 11.455 

WLON - 69 . 7b 1 

STO - SI N( COLAT/RAD) 

CTO = COS(COLAT/RAD) 

GO TO (1,2), I 

1 CTG = S I N( GLAT/RAD ) 

STG = COS(GLAT/RAD) 

CTD - CTO-CTG ♦ STO*STG*COS( (GLON+WLON)/RAD) 

STO - S()RT( 1. -CTD**2 ) 

DLAT - RAD • ATAN( CTD/STD ) 

CLD - ( CTO-CTD-CTG) / (STO-STD) 

SLGLO - SIN((GLON+WLON)/RAD) 

SLU = STG * SLGLO / STD 

DLON - 180. - RAU-ATAN2 (SLD,-CLD) 

GO TO 3 

2 CTD - S I N( DLAT/RAD ) 

STD = COS (DLAT/RAD) 

CTG - CTO*CTD - STO*STD*COS(DLUN/RAD) 

STG - SQRT( 1. -CTG**2) 

GLAT - RAD * ATAN( CTG/STG ) 

CLGLO - ( CTD-CTO-CTG) / (STO-STG) 

SLGLO * STD * SI N(D«.ON/RAD) / STG 
GLUN = RAD*ATAN2(SLGL0, CLGLO) - WLON 
I F ( GLON . LT. 0. ) GLON - GLON*3bO. 

5 COSD - (CTO-CTD-CTG) / ( STD-STG) 

SIND - STO * SLGLO / STD 
D - RAD * ATAN2( SIND, COSD) 

RETURN 

END 

SUBROUTINE Dl PFLD( TM, RKM, ST, CT, SL, CL, BR, BT, BP, B ) 

FIELD MODEL IS IGRF 1965.0 IN DIPOLE COORDINATES (MEAD, 1970). 

TM IS TIME IN YEARS (E.G. 1969.2) FOR WHICH THE FIELD IS DESIRED. 
INPUT ARE GEOCENTRIC DISTANCE IN KM (RKM) AND SINE AND COSINE OF 
DIPOLE COLATITUDE ($T,CT) AND EAST DIPOLE LONGITUDE W/R/T THE 
DIPOLE MERIDIAN PASSING THROUGH THE SOUTH GEOGRAPHIC POLE (SL,CL). 
OUTPUT BT AND BP (GAMMAS) ARE ALSO WITH RESPECT TO DIPOLE AXIS. 
DIMENSION F(9,9),FT(9,9), SHMI T( 9, 9 ) , CONST (9,9),SP(9),CP(9), 

1 P(9,J),DP(9,9),G(9,9) 

DATA SUM I T ( 1, 1),P( 1, 1),CP( 1),DP( 1,1), SP( 1) / -1. , 2*1. , 2*0. / 

DATA TZERO, NMAX /1965.0,9/ 

DATA F/0., -309 5 3., -6 18., 906., 8 37., -14 0.., 48., 69., 10., 0., 0., 2997., 
1-1238. , -490. , 90. ,-17., -40. ,-l., 2255. ,481. ,-1875. ,-1052. ,15. ,-322., 

2 60., 8., 12., -1758., 1170., -485., -546., 311., 5 3., 18 5., -14., -11., 962., 

3 17b., -39., -3 12., -317., -138., -151., 15., -14., 344., -49., 170., 103., 

4 4b., 7., 22., 3., 4., 7., -55., 7 8., -56., -29., -95., -4 8., 27., -2b., -39., 

5 bl.,-9.,36.,4.,-3.,10.,5.,4.,12.,-2.,1.,4.,7.,10.,11.,9.,-10./ 
DATA FT/0.,16.,-18.8,-3.4,-. 3, -. 5, . 3, . 6, -. 1, 7. 4, 2 . , 20. 9, -6 . 8, 1. 1, 

1 -3.9,0.,!., -.3, -1.4, 12. 2, 8. 8, 1.3, 1.2, -2., -1.3,0., -.3, -8., 4. 6, 8., 
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3 .3#1.4/ B .7/ B .4..8. B .2. B .l/.3**»4. B .3/"2*3/1.7*".2/.3f B .8/"*2*.6. 

4 B li# B tl/»l/0*#*3f**5/ a iS/i2/«2/*7|«li B tl/ B »2# B *l#t2# • 1 # . 4 . *2/ 

I F <SHMt T( 2,1) .EQ. -1. ) GO TO 3 

DO 2 N«2,NMAX 

SHMI T(N, 1) - SHMI T( N-l, 1) * <2*N-3) / (N-l) 

DO 2 M«1,N 
1 F CM- EQ. 1 ) GO TO 1 
JJ » 1 

I F (M. EQ. 2 ) JJ « 2 

SHMIT(N,M) - SHMIT(N,M-1) * SQRT( FLOAT ( (N-M+l)*JJ)/(N+M-2) ) 
F(M-1,N) • F(M-1,N) * SHMI T( N,M) 

FT(M-1,N) ■ FT(M-1, N) * SHMI T( N,M) 

1 F(N,M) - F(N,M) * SHMIT(M,M) 

FT(N,M) « FT(N,M) * SHMI T(N,M) 

2 CONST( N,M) - FLOAT ( (N-2)**2-(M-l)**2) / < ( 2*N-3 )*< 2*N- 5 ) ) 

3 IF(TM.EQ.TMLAST) GO TO 5 
TMLAST ■ TM 

T - TM-TZERO 
DO 4 N-2,NMAX 
DO 4 M»1,N 

G(N,M) « F(N,M) ♦ T*FT(N,M) 

IF(M.EQ.l) GO TO 4 

G(M-1,N) - F(M-1,N) ♦ T*FT(M-1, N) 

4 CONTINUE 

C CALCULATION FOR EACH POINT USUALLY BEGINS HERE. 

5 CONTINUE 
SP(2) ■ SL 
CP( 2 ) - CL 

DO 6 M«3, NMAX 

S P ( M ) • SP(2)*CP(M-1) ♦ CP( 2 )*SP(M-1 ) 

6 CP(M) - CP( 2)*CP(M-1) - SP(2)*SP(M-1) 

AOR - 6371.16 / RKM 

AR - AOR**2 
BR - 0. 

BT - 0. 

BP » 0. 

DO 11 N«2,NMAX 
AR - AOR*AR 
DO 11 M*l, N 
IF(M.EQ.N) GO TO 7 

P( N,M) - CT*P(N-1,M) - CONST(N,M)*P<N-2,M) 

DP(N/M) - CT*OP(N-l,M) - ST*P(N-1,M) - CONST(N,M)*DP(N-2,M) 

GO TO 8 

7 PCN,N) - ST*P(N-1,N-1) 

DP(N,N) ■ ST*DP< N-l, N-l ) * CT*P(N-1,N-1) 

8 PAR - P(N/M)*AR 

IF (M.EQ.l) GO TO 9 

TEMP - G(N,M)*CP<M) ♦ G(M-1, N)*SP(M) 

BP - BP - (M-l) * PAR * (G(N,M)*SP(M)-G(M-1,N)*CP(M)) 

GO TO 10 

9 TEMP - G(N,M) 

10 BT - BT ♦ AR*OP(N,M)*TEMP 

11 BR - BR - PAR*N*TEMP 
BP - BP/ST 

B « SQRT(BR**2*8T**2*BP**2) 

RETURN 

END 
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